// Copyright 2022 Jeroen van Nugteren

// Permission is hereby granted, free of charge, to any person obtaining a copy of this software
// and associated documentation files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use, copy, modify, merge, publish,
// distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:

// The above copyright notice and this permission notice shall be included in all copies or
// substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
// OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
// CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

// include header
#include "dfplasma.hh"

#include "rat/common/error.hh"
#include "dfpolygon.hh"

// code specific to Rat
namespace rat{namespace dm{

	// constructor
	DFPlasma::DFPlasma(){
		set_name("Plasma");
	}

	// constructor with shape input
	DFPlasma::DFPlasma(
		const rat::fltp a,
		const rat::fltp r0, 
		const rat::fltp delta, 
		const rat::fltp kappa, 
		const bool center_r0){
		
		set_name("plasma");
		set_a(a); set_r0(r0); set_delta(delta); 
		set_kappa(kappa); set_center_r0(center_r0);
	}

	// factory
	ShDFPlasmaPr DFPlasma::create(){
		return std::make_shared<DFPlasma>();
	}

	// factory
	ShDFPlasmaPr DFPlasma::create(
		const rat::fltp a, 
		const rat::fltp r0, 
		const rat::fltp delta, 
		const rat::fltp kappa, 
		const bool center_r0){
		return std::make_shared<DFPlasma>(a,r0,delta,kappa,center_r0);
	}


	// set plasma parameters
	void DFPlasma::set_r0(const fltp r0){
		r0_ = r0;
	}

	void DFPlasma::set_center_r0(const bool center_r0){
		center_r0_ = center_r0;
	}

	void DFPlasma::set_a(const fltp a){
		a_ = a;
	}

	void DFPlasma::set_delta(const fltp delta){
		delta_ = delta;
	}

	void DFPlasma::set_kappa(const fltp kappa){
		kappa_ = kappa;
	}




	// set plasma parameters
	fltp DFPlasma::get_r0() const{
		return r0_;
	}

	bool DFPlasma::get_center_r0() const{
		return center_r0_;
	}

	fltp DFPlasma::get_a() const{
		return a_;
	}

	fltp DFPlasma::get_delta() const{
		return delta_;
	}

	fltp DFPlasma::get_kappa() const{
		return kappa_;
	}


	// get bounding box
	arma::Mat<fltp>::fixed<2,2> DFPlasma::get_bounding() const{
		// allocate
		arma::Mat<fltp>::fixed<2,2> Mb;
		
		// bounding in x
		Mb.col(0) = arma::Col<fltp>::fixed<2>{-a_,a_};

		// bounding in y
		Mb.col(1) = arma::Col<fltp>::fixed<2>{-kappa_*a_, kappa_*a_};

		// center
		if(!center_r0_)Mb.col(0) += r0_;

		// return box
		return Mb;
	}

	// perimeter function
	ShPerimeterPr DFPlasma::create_perimeter(const fltp delem) const{
		// calculate number of elements
		const arma::uword num_elements = std::max(arma::uword(std::ceil(2*arma::Datum<fltp>::pi*a_/delem)),2llu);

		// create node coordinates on a circle
		const arma::Row<fltp> theta = arma::linspace<arma::Row<fltp> >(
			0,((static_cast<fltp>(num_elements)-1.0)/num_elements)*arma::Datum<fltp>::tau,num_elements);
		const arma::Mat<fltp> Rn = arma::join_vert(
			a_*arma::cos(theta + delta_*arma::sin(theta)),
			kappa_*a_*arma::sin(theta));

		// create element matrix
		const arma::Mat<arma::uword> n=arma::join_vert(
			arma::regspace<arma::Row<arma::uword> >(0,Rn.n_cols-1),
			arma::shift(arma::regspace<arma::Row<arma::uword> >(0,Rn.n_cols-1),-1,1));
		
	 	// create perimeter
		return Perimeter::create(Rn,n);
	}

	// distance function
	arma::Col<fltp> DFPlasma::calc_distance(const arma::Mat<fltp> &p) const{
		// convert to carthesian coordinates
		const arma::uword num_points = 180;
		const arma::Col<fltp> theta = arma::linspace<arma::Col<fltp> >(0,2*arma::Datum<fltp>::pi,num_points);
		
		// calculate v
		arma::Mat<fltp> v = arma::join_horiz(
			a_*arma::cos(theta + delta_*arma::sin(theta)),
			kappa_*a_*arma::sin(theta));
		
		// center
		if(!center_r0_)v.col(0) += r0_;

		// number of points and number of sections
		arma::uword np = p.n_rows;
		arma::uword nv = v.n_rows;

		// distance matrix
		arma::Mat<fltp> ds(np,nv);

		// walk over sections
		for(arma::uword i=0,j=nv-1;i<nv;j=i++){
			// get vertices
			const arma::Row<fltp>::fixed<2> v1 = v.row(i);
			const arma::Row<fltp>::fixed<2> v2 = v.row(j);
				
			// calculate edge vector
			const arma::Row<fltp>::fixed<2> dv = v2-v1;

			// walk over points
			for(arma::uword ip=0;ip<np;ip++){
				// get point
				const arma::Row<fltp>::fixed<2> myp = p.row(ip);

				// calculate vector from first vertex to point
				const arma::Row<fltp>::fixed<2> dp = myp - v1;

				// calculate dot products
				const fltp c1 = arma::accu(dv%dp);
				const fltp c2 = arma::accu(dv%dv);
				
				// distance to projected point on line
				ds(ip,i) = std::sqrt(arma::accu(arma::square(
					myp-(v1+std::min(RAT_CONST(1.0),std::max(RAT_CONST(0.0),c1/c2))*dv))));
			}
		}

		// check which points are inside the polygon
		const arma::Col<fltp> sign = 2*arma::conv_to<arma::Col<fltp> >::from(DFPolygon::inpolygon(v,p)) - 1;

		// check output
		assert(sign.is_finite());

		// signed distance
		arma::Col<fltp> sds = -arma::min(ds,1)%sign;

		// negative distance inside
		return sds;
	}

	// validity check
	bool DFPlasma::is_valid(const bool enable_throws) const{
		// check input again
		if(r0_<=0){if(enable_throws){rat_throw_line("major radius must be larger than zero");} return false;}
		if(a_<=0){if(enable_throws){rat_throw_line("minor radius must be larger than zero");} return false;}
		// if(delta_<=0){if(enable_throws){rat_throw_line("triangularity must be larger than zero");} return false;}
		if(kappa_<=0){if(enable_throws){rat_throw_line("kappa must be larger than zero");} return false;}
		return true;
	}

	// get type
	std::string DFPlasma::get_type(){
		return "rat::dm::dfplasma";
	}

	// method for serialization into json
	void DFPlasma::serialize(
		Json::Value &js, rat::cmn::SList &list) const{
		
		// parent
		DistFun::serialize(js,list);

		// type
		js["type"] = get_type();

		// properties
		js["center_r0"] = center_r0_;
		js["ro"] = r0_; 
		js["a"] = a_; 
		js["delta"] = delta_; 
		js["kappa"] = kappa_;
	}

	// method for deserialisation from json
	void DFPlasma::deserialize(
		const Json::Value &js, rat::cmn::DSList &list, 
		const rat::cmn::NodeFactoryMap &factory_list, 
		const boost::filesystem::path &pth){

		// parent
		DistFun::deserialize(js,list,factory_list,pth);

		// properties
		r0_ = js["ro"].ASFLTP(); 
		a_ = js["a"].ASFLTP();
		delta_ = js["delta"].ASFLTP(); 
		kappa_ = js["kappa"].ASFLTP();
		center_r0_ = js["center_r0"].asBool();
	}

}}